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In  graphene,  as  in  most  metals,  electron-electron  interactions  renormalize  the  properties  of  electrons  but 
leave  them  behaving  like  noninteracting  quasiparticles.  Many  measurements  probe  the  renormalized 
properties  of  electrons  right  at  the  Fermi  energy.  Uniquely  for  graphene,  the  accessibility  of  the  electrons 
at  the  surface  offers  the  opportunity  to  use  scanned  probe  techniques  to  examine  the  effect  of  interactions 
at  energies  away  from  the  Fermi  energy,  over  a  broad  range  of  densities,  and  on  a  local  scale.  Using 
scanning  tunneling  spectroscopy,  we  show  that  electron  interactions  leave  the  graphene  energy  dispersion 
linear  as  a  function  of  excitation  energy  for  energies  within  ±200  meV  of  the  Fermi  energy.  However,  the 
measured  dispersion  velocity  depends  on  density  and  increases  strongly  as  the  density  approaches  zero 
near  the  charge  neutrality  point,  revealing  a  squeezing  of  the  Dirac  cone  due  to  interactions. 


DOI:  10. 1103/PhysRevLett.  109. 116802 


PACS  numbers:  73.22.Pr,  68.37.Ef,  73.20.At 


In  the  absence  of  interactions  the  low  energy  excitations 
in  graphene  are  described  by  massless  Dirac  quasiparticles 
with  a  linear  energy-momentum  dispersion  E  =  hvk , 
where  the  constant  of  proportionality  is  the  carrier  group 
velocity  v.  When  many-body  interactions  are  included,  the 
energy  dispersion  can  change  [1-12].  In  general,  this 
change  depends  on  the  parameters  [n,  k,  r  J,  where  n  is 
the  carrier  density,  k  is  the  momentum  relative  to  the 
Brillouin  zone  corner,  and  rs  is  the  interaction  parameter 
or  coupling  constant  describing  the  relative  strength  of  the 
Coulomb  interactions  and  is  given  by  the  ratio  of  potential 
to  kinetic  energies  [13].  In  graphene,  both  the  kinetic  and 
potential  energies  scale  as  the  square  root  of  the  density,  so 
that  rs  is  independent  of  density.  Therefore,  in  graphene, 
the  dependence  of  the  dispersion  renormalization  on  rs  and 
n  can  be  separated.  For  monolayer  graphene  rs  is  given  by 
rs  =  e2/47Ts0smhv ,  where  e  is  the  electron  charge,  h  is 
Planck’s  constant  divided  by  2i 7,  e0  is  the  permittivity  of 
free  space,  and  sm  is  the  effective  dielectric  constant  of  the 
medium  in  which  the  graphene  sheet  is  embedded. 

Previously,  measurements  of  the  graphene  dispersion  re¬ 
normalization  have  either  examined  solely  the  density  de¬ 
pendence  [1 1]  or  the  dependence  on  energy  [6,9,10]  but  not 
both  over  wide  ranges  of  density  and  energy  due  to  limita¬ 
tions  in  experimental  techniques.  Transport  measurements, 
which  are  sensitive  to  the  behavior  of  electrons  at  the  Fermi 
energy  ( E¥ ),  have  shown  the  dispersion  velocity  at  EF  in¬ 
creases  at  low  density  [11].  Angle-resolved  photoemission 
spectroscopy  has  examined  the  energy  dependence  of  the 


dispersion  and  found  deviations  from  linearity  at  high  ener¬ 
gies  below  approximately  —0.4  eV  [9].  Angle-resolved  pho¬ 
toemission  spectroscopy  measurements  are  complicated  by 
the  appearance  of  plasmonic  and  phonon  structures  which 
can  distort  the  spectrum  at  the  Dirac  point  [6],  and  they  have 
not  been  able  to  examine  the  density  dependence  in  detail  due 
to  the  lack  of  large  area  gated  graphene  devices.  Scanning 
tunneling  spectroscopy  measurements  of  back-gated  devices 
allows  the  graphene  dispersion  to  be  examined  as  a  function 
of  energy  away  from  the  Fermi  energy  [8,14],  so  that  the 
renormalization  can  be  determined  separately  as  a  function 
of  excitation  energy  and  density. 

For  excitations  at  the  Fermi  level,  i.e.,  at  k  =  kF, 
many-body  theory  predicts  that  the  renormalized  spectrum 
of  graphene  can  be  characterized  by  a  velocity  u*  given 
by  [1-3] 


v  tt 


(1) 


The  Fermi  velocity  enhancement  described  in  Eq.  (1)  was 
calculated  within  the  random  phase  approximation,  where 
v  is  the  bare  dispersion  velocity  in  the  absence  of  inter¬ 
actions  and  nc  is  a  density  corresponding  to  the  ultraviolet 
cutoff  energy,  which  is  ~  3  eV  [2,11].  The  2nd  and  3rd 
terms  in  Eq.  (1)  result  from  intraband  and  interband  con¬ 
tributions,  respectively.  The  expression  in  Eq.  (1)  describes 
the  reshaping  of  the  Dirac  cone  at  energies  equal  to  the 
Fermi  energy  EF.  Velocity  enhancements  obtained  from 
Shubnikov-de  Haas  oscillations  in  transport  measurements 
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have  considered  only  the  interband  contribution  in  Eq.  (1) 
and  found  the  need  for  an  additional  density-dependent 
parameter  to  get  agreement  with  experiment  [11].  The 
question  on  how  the  dispersion  velocity  depends  on  energy 
over  a  broad  range  of  densities  remains  open. 

In  this  Letter,  we  present  the  first  experimental  measure¬ 
ments  of  the  Landau  level  (LL)  spectroscopy  of  graphene  on 
hexagonal  boron  nitride  (hBN)  spacer  layers.  The  lower 
disorder  in  graphene  with  hBN  spacer  layers  allows  a 
significant  improvement  in  LL  lifetimes,  with  the  appear¬ 
ance  of  many  LLs  over  a  wide  energy  range.  This  enhanced 
spectrum,  compared  to  previous  results  without  hBN  [8,14], 
allows  us  to  separate  the  energy  and  density  contributions  to 
the  renormalization  of  the  graphene  band  structure.  We 
show  how  to  quantitatively  examine  the  Landau  level  tun¬ 
neling  spectrum  taking  into  account  probe  tip  and  graphene 
quantum  capacitances  and  to  extract  interaction  driven 
velocity  renormalization  with  good  accuracy.  We  find  that 
interactions  do  not  significantly  distort  the  Dirac  cones  at 
low  energies.  Instead,  they  preserve  the  linear  dispersion 
while  parametrically  renormalizing  the  dispersion  velocity 
at  a  given  fixed  density.  The  measured  renormalized  veloc¬ 
ity  is  satisfactorily  described  by  the  random  phase  approxi¬ 
mation  theory  incorporating  the  interaction  strength  rs  and 
the  ultraviolet  cutoff  density  nc  [2,12]  described  by  Eq.  (1) 
without  the  need  of  additional  parameters. 

The  experiments  were  performed  in  a  custom  designed 
cryogenic  scanning  tunneling  microscope  system  operat¬ 
ing  at  4  K  with  applied  magnetic  fields  from  0  to  8  T  [15]. 
The  graphene  devices  were  fabricated  by  the  method  de¬ 
tailed  in  Dean  et  al.  [16].  A  heavily  doped  silicon  substrate 
was  used  as  a  back  gate  to  control  the  carrier  density  n  in 
the  sample  by  applying  back  gate  voltage  VG.  In  our 
experiments,  a  bias  voltage  VB  is  applied  to  the  graphene, 
and  the  tunneling  current  I  is  measured  from  the  tip  to 
obtain  tunneling  spectra  dl/dV  as  a  function  of  VB,  VG, 
and  applied  magnetic  field  B.  Topography  measurements 
of  graphene  on  hBN  [Fig.  1(a)]  are  characterized  by  a 
moire  pattern  formed  by  the  relative  rotation  of  the  gra¬ 
phene  sheet  with  respect  to  the  underlying  hBN  crystal  as 
observed  in  previous  studies  [17,18].  The  topographic 
modulation  in  Fig.  1(a)  is  consistent  with  the  thin  BN 
spacer  layer  (5  nm)  used  in  these  studies  [16]. 

Upon  applying  a  magnetic  field,  the  tunneling  spectra 
develops  sharp  peaks  as  a  function  of  VB  as  shown  in 
Fig.  1(b),  indicating  the  formation  of  the  LLs  [14,19-21]. 
The  reduced  disorder  in  graphene  on  hBN  is  evident  by  the 
high  fidelity  of  the  LL  formation  compared  to  previous 
measurements  on  Si02  [8,14].  The  LL  index  N  can  be 
identified  from  the  peak  distance  in  bias  VB  compared  to 
the  expected  LL  energy  spectra  proportional  to  |7V| 1//2.  The 
widths  of  the  LLs,  determined  by  the  lifetime  of  the  Dirac 
quasiparticles,  approaches  those  measured  in  the  low- 
disorder  graphene  grown  epitaxially  on  SiC  [19,21]. 
Figure  1(c)  shows  the  map  of  the  N  =  —  1  LL  peak 
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FIG.  1  (color  online).  Local  probing  of  graphene  on  hBN. 

(a)  STM  topographic  image  of  graphene  on  hBN.  The  moire 
pattern  with  a  unit  cell  length  of  4.5  nm  corresponds  to  a  rotation 
of  3.1°  of  the  graphene  lattice  relative  to  the  hBN  crystal. 

(b)  Differential  conductance  spectra  of  graphene  on  hBN  at 
VG  =  -30  V.  (c)  Spatial  variation  of  the  N  =  -1  LL  peak 
position  measured  in  the  same  area  as  (a),  reflecting  the  local 
disorder  potential  variation.  The  positions  A  and  B  (white  dots) 
denote  electron  and  hole  puddle  locations  where  measurements 
are  reported  in  Figs.  2-4. 


position,  which  reflects  the  spatial  variation  of  disorder 
potential.  Within  the  area  of  interest,  two  local  potential 
extrema  are  indicated  as  A  and  B.  At  low  density,  A  becomes 
an  electron  puddle  (blue)  while  B  a  hole  puddle  (yellow). 

The  tunneling  spectra  in  Fig.  1(b)  reflect  the  graphene 
LL  density  of  states  at  a  gate  voltage  of  VG  =  —30  V, 
which  induces  additional  holes  to  lower  the  Fermi  level  by 
~  200  meV  with  respect  to  the  Dirac  point.  Repeating  the 
same  tunneling  spectra  measurements  while  varying  the 
back  gate  voltage,  thus  controlling  the  relative  positions  of 
the  Fermi  level  and  the  Dirac  point,  we  can  obtain  a 
complete  data  set,  which  we  refer  to  as  a  “gate  map.” 
The  resulting  tunneling  spectra,  represented  in  a  two- 
dimensional  plot  of  VB  and  VG,  are  shown  in  Figs.  2(a) 
and  2(b),  for  2  and  5  T,  respectively.  The  LL  spectra 
variation  as  a  function  of  energy  and  density  can  be  inves¬ 
tigated  by  tracing  peaks  of  the  tunneling  spectra  (bright 
curve  traces)  in  these  gate  maps.  At  a  lower  magnetic  field 
[Fig.  2(a)],  the  smooth  curvature  of  the  LLs  with  gate 
voltage  approximately  corresponds  to  the  energy  variation 
of  the  charge  neutrality  point  (Dirac  point),  which  varies  as 
Ed  oc  n  V2.  In  larger  fields  [Fig.  2(b)],  LL  formation  leads 
to  the  familiar  stair  case  pattern  in  the  gate  maps,  previ¬ 
ously  seen  in  graphene  on  Si02  [8,14],  or  GaAs  two- 
dimensional  electron  gases  [22].  Every  LL  becomes  pinned 
at  the  Fermi  level  until  it  fills  with  carriers,  and  upon  filling 
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the  next  unfilled  LL  makes  a  quick  transition  to  become 
pinned  at  EF. 

To  accurately  analyze  the  energy  dispersion,  we 
simulate  the  single  particle  properties  in  the  gate  maps 
[Fig.  2(c)]  by  calculating  the  tunneling  conductance  spec¬ 
tral  map  using  a  capacitor  model  [23]  that  includes  the 
graphene-back-gate  capacitance  (CG),  graphene-probe  ca¬ 
pacitance  (CP),  the  graphene  quantum  capacitance,  and  a 
constant  velocity  v.  The  simulation  in  Fig.  2(c)  shows  that 
the  LL  staircase  transitions  do  not  occur  at  fixed  VG  but  at 
different  values  as  VB  changes,  appearing  as  a  line  with  an 
angle  in  the  gate  map.  This  angle  is  due  to  the  local  gating 
of  the  graphene  by  the  bias  potential  between  the  probe  tip 
and  graphene  and  is  determined  by  the  ratio  of  the  back 
gate  and  probe  capacitances,  CG/CP.  The  constant  density 
(or  constant  chemical  potential)  axis  is  therefore  along  the 
LL  transitions  given  by  this  capacitance  ratio,  as  indicated 
by  the  tilted  red  line  in  Fig.  2(b).  For  a  better  comparison, 
we  simulate  the  LL  spectra  in  the  VG  and  VB  maps  by 
considering  a  fixed  v  and  a  smoothly  varying  CG/CP  ratio 
determined  from  the  experimental  gate  maps  [24],  as 
shown  in  Fig.  2(c).  The  simulation  result  matches  well 
the  experimental  LL  transitions  and  indicates  that  the 
simulation  captures  the  single  particle  physics.  Departure 
from  single  particle  behavior  due  to  many-body  effects  is 
seen  by  a  comparing  the  simulation  with  the  experiment  at 
low  densities.  Deviations  from  the  constant- velocity  simu¬ 
lation  are  illustrated  by  overlapping  the  simulation  from 
Fig.  2(c)  onto  Fig.  2(b)  (yellow  overlay)  in  the  regions  of 
high  density  and  near  the  charge  neutrality  point.  Very 
good  agreement  is  seen  in  matching  the  LL  energy  peak 
positions  at  high  density.  However,  the  simulated  LL  peak 
positions  underestimate  the  experimental  positions  and 
significantly  deviate  at  low  density  indicating  that  the 
dispersion  velocity  is  larger  close  to  the  charge  neutrality 
point.  For  example,  the  simulated  peak  position  of  LL2  lies 
in  the  minimum  of  the  spectrum  between  the  peaks  of  LLi 
and  LL2  at  zero  gate  voltage  [see  the  white  arrow 
in  Fig.  2(b)]. 

To  accurately  determine  the  interaction  driven  velocity 
enhancements,  we  define  a  constant  density  (or  chemical 
potential)  axis  in  the  (VB,  VG)  space  using  the  transitions  of 
the  LLs,  which  is  determined  from  the  CG/CP  ratio  [24]. 
Therefore  the  LL  energy  EN  can  be  obtained  at  a  fixed 
density  for  the  Vth  LL  (LL#)  by  measuring  the  LL  spectral 
peak  positions  in  the  plateau  regions  along  lines  of 
constant  density  [the  tilted  red  line  in  Fig.  2(b)]  thereby 
converting  dI/dV(VB,  VG)  to  dI/dV(E,  n).  A  comparison 
of  the  dl/dV  spectra  at  constant  VG  vs  constant  density  n 
[Fig.  3(a)]  shows  that  the  spectrum  at  constant  gate  voltage 
overestimates  the  energy  scale  by  as  much  as  5%  leading  to 
a  corresponding  error  in  an  estimate  of  the  dispersion 
velocity.  The  importance  of  examining  the  correct  LL# 
peak  energies  at  constant  density  is  compounded  by  the 
fact  that  the  energy  scale  error  at  constant  gate  voltage 
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FIG.  2  (color  online).  Gate  mapping  tunneling  spectroscopy  of 
the  Landau  level  density  of  states  of  graphene  on  hBN.  Each  map 
is  built  of  individual  dl/ dV  vs  VB  spectra  taken  at  multiple  fixed 
VG.  The  color  scale  is  the  dl/dV  magnitude,  (a)  dl/dV  gate  map 
at  B  =  2  T  obtained  in  the  hole  puddle  labeled  B  in  Fig.  1(c). 
(b)  dl/dV  gate  map  spectra  at  B  =  5  T.  The  tilted  red  line 
indicates  an  axis  of  constant  density  (chemical  potential)  follow¬ 
ing  the  transitions  of  the  LLs  as  discussed  in  the  text,  while  the 
vertical  blue  line  indicates  constant  VG.  The  yellow  overlays  are 
portions  of  the  simulation  in  (c).  (c)  Simulation  of  the  dl/dV 
gate  map  spectra  at  B  =  5  T  with  a  constant  dispersion  velocity 
of  1.1  X  106  m/s  and  the  probe-sample  capacitance  determined 
from  the  slope  of  the  LL  transitions  in  (b)  [24].  Portions  of  the 
simulation  are  overlaid  in  (b)  (yellow  lines)  showing  agreement 
at  high  densities  and  deviations  at  low  densities  from  using  a 
constant  velocity.  The  white  arrow  points  to  the  simulated  LL2 
position  which  lies  in  between  the  experimental  LLj  and  LL2 
peaks  at  low  density. 
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FIG.  3  (color  online).  Landau  level  dispersion,  (a)  dl/dV 
spectra  (red  line)  at  constant  density  obtained  along  the  red 
line  in  Fig.  2(b).  A  fit  of  LL  peaks  N  =  1  to  N  =  6  is  shown 
in  green  using  a  series  of  Lorentzians.  (Inset)  A  comparison  of 
the  dll dV  spectra  at  constant  gate  voltage  (VG  =  12  V)  (blue 
line)  and  constant  density  (red  line),  obtained  along  the  lines 
indicated  in  Fig.  2(b).  The  spectra  at  constant  gate  voltage 
overestimates  the  energy  scale  by  5%  at  this  density, 
(b)  Determination  of  the  graphene  dispersion  velocity  from  the 
LL  peak  energy  positions  (obtained  along  lines  of  constant 
density  in  the  gate  maps)  for  the  electron  puddle  A  at  B  = 
2  T.  The  peak  positions  are  plotted  as  a  function  of  the  square 
root  of  the  Landau  orbital  index  N  and  magnetic  field  B.  A  linear 
fit  (solid  lines)  is  used  to  determine  the  dispersion  velocity. 
(Inset)  Residuals  from  the  linear  fit  showing  very  good  linearity 
in  the  LL  dispersion. 

varies  with  density,  since  the  CG/CP  ratio  varies.  This  will 
result  in  errors  as  large  30%  in  determining  rs  parameters 
from  the  velocity  incorrectly  determined  at  constant  gate 
voltage.  Below,  we  determine  the  many-body  corrections 
to  the  renormalized  velocity  by  correctly  analyzing  the  LL 
spectra  at  constant  density  with  a  correct  energy  scale  [24]. 

In  graphene,  the  linear  dispersion  yields  the  LL  spectra: 
En  =  Ed  +  sgn(N)s/2ehv2NB.  The  LL  spectroscopy  pre¬ 
sented  in  Fig.  2  thus  allows  us  to  probe  the  dispersion 
velocity  at  fixed  density  and  magnetic  field  by  inspecting 
En  for  different  N  [Fig.  3(a)].  Figure  3(b)  shows  EN-ED  vs 
^JNB  for  different  density  values.  Here  ED  is  obtained  from 
En= o  at  a  given  density.  A  remarkable  feature  in  Fig.  3(b) 
is  that  the  LL  dispersion  is  highly  linear  [see  the  inset  in 
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FIG.  4  (color  online).  Renormalized  graphene  dispersion, 
(a)  Schematic  of  the  Dirac  cone  variation  as  a  function  of 
density.  The  data  presented  in  this  Letter  show  that  the  graphene 
energy-momentum  dispersion  remains  linear  at  low  energy  in  the 
presence  of  electron  interactions,  while  the  Dirac  cone  angle, 
which  is  inversely  proportional  to  the  velocity,  decreases  (gets 
squeezed)  at  low  density,  (b)  Renormalized  velocity  determined 
from  the  linear  fitting  of  LL  peak  positions  as  described  in 
Fig.  3(b)  as  a  function  of  density  for  the  electron  (red  symbols) 
and  hole  (blue  symbols)  puddle  locations  A  and  B  at  B  =  2  T. 
The  solid  line  shows  a  fit  to  the  combined  data  from  puddle  A 
and  B  using  Eq.  (1)  with  v  =  (0.957  ±  0.003)  X  106  m/s  and 
rs  =  (0.69  ±  0.03)  [25]. 


Fig.  3(b)]  at  a  fixed  density  but  changes  slope  as  the  density 
decreases.  The  energy  dispersion  is  thus  described  by  the 
Dirac  cone  in  this  energy  range,  but  the  cone  is  squeezed  at 
low  density,  as  schematically  illustrated  in  Fig.  4(a). 

The  velocity  extracted  from  the  LL  dispersion  fits  is 
shown  in  Fig.  4(b)  as  a  function  of  density  for  the  electron 
and  hole  puddle  locations,  A  and  B.  We  observe  a  signifi¬ 
cant  increase  in  velocity  near  the  charge  neutrality  point. 
We  fit  the  renormalized  velocity  using  Eq.  (1).  The  solid 
line  in  Fig.  4(b)  shows  a  two-parameter  fit,  to  the  combined 
data  of  puddles  A  and  B ,  with  a  bare  velocity  v  =  (0.957  ± 
0.003)  X  106  m/s  and  an  interaction  parameter  rs  = 
0.69  ±  0.03  [25].  The  bare  velocity  is  in  agreement  with 
previous  measurements  by  a  large  number  of  different 
experiments  [19,21,26-29].  The  fit  parameter  rs  is  equiva¬ 
lent  to  a  measurement  of  the  effective  dielectric  constant  of 
£BN/sio2  =  5.3  ±  0.3  [25],  using  sm  =  (eBN/sio2  +  l)/2. 

In  summary,  we  present  measurements  of  the  LL  density 
of  states  of  graphene  on  hBN  as  a  function  of  spatial 
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position,  density,  and  magnetic  field.  Tunneling  spectros¬ 
copy  at  energies  away  from  E¥  allow  us  to  separately 
analyze  the  graphene  dispersion  dependence  on  excitation 
energy  and  density.  The  dispersion  of  the  LLs  as  a  function 
of  orbital  index  shows  that  the  linear  graphene  energy- 
momentum  dispersion  is  retained  at  low  energies  even  as 
the  charge  neutrality  point  is  approached.  However,  electron 
interactions  cause  an  increase  of  the  dispersion  velocity, 
effectively  squeezing  the  Dirac  cone,  with  decreasing  den¬ 
sity  in  agreement  with  recent  theoretical  predictions  [12]. 
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